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Nomenclature 


A 

A* 

CL 

C p 

Cv 

err 

E 

F 

9c 

ic 

Icmd 

K 

K p 

Ki 

Mwt 

M a 

N 

P 

P cmd 
PI 
Q 
R 

strk 

t 

T 

tau 

V 
w 
W 
x 
X 

Y 

Ycmd 

Ydot 


area 

minimum flow area 

effective cylinder length 

molar average constant pressure specific heat 

molar average constant volume specific heat 

error between process and set point 

energy 

molar flow 

conversion constant 

initial condition 

integral command 

shock loss coefficient 

controller proportional gain 

controller integral gain 

molecular weight 

mach number 

moles 

pressure 

proportional command 
proportional integral 
heat transfer 
universal gas constant 
piston stroke 
time 

temperature 
valve time constant 
volume 

compressor angular speed 
work 

displacement 

concentration heavy gas in mixture 
valve position 
command 
valve velocity 


in 
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P 

0 


molar average specific heat ratio 
density 

phase angle between cylinders 


subscripts 


air 

air 

cond 

condensed molar flow 

cv 

control volume 

heat 

heat transferred 

hg 

heavy gas 

in 

inlet conditions 

out 

outlet conditions 

set 

set point 

t 

stagnation conditions 

1.2 

states 1 and 2 


iv 



Abstract 


A mathematical model of a multi-stage compressor with variable molecular 
weight flow medium is derived. The modeled system consists of a five stage, six 
cylinder, double acting, piston type compressor. Each stage is followed by a 
water cooled heat exchanger which serves to transfer the heat of compression 
from the gas. A high molecular weight gas (CFC-12) mixed with air in varying 
proportions is introduced to the suction of the compressor. Condensation of the 
heavy gas may occur in the upper stage heat exchangers. The state equations 
for the system are integrated using the Advanced Continuous Simulation 
Language (ACSL) for determining the system's dynamic and steady state 
characteristics under varying operating conditions. 

1.0 Introduction 

The Transonic Dynamics Tunnel (TDT) is a large scale wind tunnel located at 
NASA Langley Research Center. This tunnel utilizes a high molecular weight 
gas for improved testing characteristics such as higher Reynolds number, lower 
speed of sound as well as lower power requirements. Currently the tunnel uses 
CFC-12 (R12) as the test medium. A major project is under way to convert the 
heavy gas medium from R12, production of which will cease December 1995, to 
HFC-134a (R134a). However, it was recognized that R134a is combustible 
under certain conditions of pressure, temperature, and concentration with air. 1 

During down times or model changes, the tunnel and plenum medium may 
consist of pure air, pure R12, or any ratio of the two gases. The R12 is 
reclaimed by compressing and cooling the mixture which condenses the heavy 
gas for storage and future use. During the process of compression, the 
temperature of the gas mixture increases. Situations may occur during the 
operation of the compressor which result in discharge temperatures exceeding 
500°F, the lower limit auto-ignition temperature for R134a/air mixtures. 1 

A control strategy has been devised to limit the compression ratio across any 
stage of the compressor, thereby limiting the maximum discharge temperature. 

1 Babcock, D. A., Bruce, R. A., "Combustibility Tests of 1,1,1,2-tetrafluoroethane in a 
Simulated Compressor Cylinder." Unpublished Report, 1995. 
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In order to evaluate the control strategy and compressor operation and 
performance, a mathematical model of the system was developed. The system 
was modeled using the Advanced Continuous Simulation Language (ACSL) 
(Reference 1) computer code to solve the system of nonlinear time dependent 
ordinary differential equations. The model utilizes a lumped volume approach, 
solving the continuity of mass and energy equations. Friction loss and 
restrictions are represented by flow resistance coefficients. Spatial variations 
are neglected which is a suitable assumption for low velocity systems. 

ACSL is a FORTRAN based simulation language convenient for modeling 
dynamic systems. Variables may be changed at run time to evaluate their 
impact on system dynamics. A choice of integration routines is available. A 
fourth order Runga-Kutta routine with fixed time step is currently used. The time 
step is based on the lowest characteristic time constant of the system. 

In addition to control strategy and compressor operation evaluation, ignition of 
the R134a/air mixture in the compressor cylinders was simulated. This 
simulation included the ignition of the gaseous mixture due to high cylinder 
pressure ratios. Relief valve dynamics and system interaction were evaluated 
to assess resultant pressures. 

Model dynamic and steady state accuracy were evaluated by comparing 
compressor data with simulations. The model was then modified for simulating 
R134a use and other modifications to the system. 

2.0 Model Description 

Figure 1 is a schematic of the modeled system. The system consists of four 
control valves, 7 fixed volumes, 12 volumes with moving boundaries and 32 
branch flows. A gaseous mixture flows into the suction volume where the 
pressure is regulated by a control valve. Flow enters the first stage of the 
machine which consists of two, double acting cylinders. Each cylinder volume 
has a suction valve which acts essentially as a check valve limiting flow to one 
direction. The gas is then compressed and discharged through discharge 
valves in the cylinder to a common manifold. The gas then enters a water 
cooled heat exchanger which removes the heat added during the compression 
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cycle. The remaining stages of the machine operate similarly to the first stage 
but are single cylinders. Typically the pressure in the third stage heat 
exchanger is high enough to condense heavy gas. A condensation line is 
shown leaving the third stage heat exchanger. In order to keep high pressure in 
the third stage, gas is recirculated from the fifth stage discharge to the third 
stage discharge via a pressure control bypass valve. Flow proceeds to the 
fourth and fifth stages of the machine where condensation may also occur. 
Discharge pressure is maintained by a PI pressure control valve at the fifth 
stage discharge. An overall bypass valve may unload the machine by routing 
discharge flow to machine inlet. This reduces overall machine pressure ratio 
and individual stage pressure ratios such that stage discharge temperatures 
remain in a safe regime. 

2.1 Volume Analysis 


A schematic of a system with moving boundaries is shown in Figure 2. On its 
downward stroke, the piston pulls a charge of gaseous mixture through the 
suction valve into the control volume. The piston then compresses and 
discharges the gas through the discharge valve into the piping system. 
Assuming adiabatic walls, the conservation equations of mass and energy are 
written and manipulated such that the system state variables are pressure and 
temperature. 


Piston position for any time is given by: 


x = 


r strk A 


sin (wt + <(>) 


Instantaneous volume is given by: 


V cv = (CL - x) 


CL- 


" strk^ 

2 J 


\ 

sin (wt + <j)) 

/ 


A 

cv 


With cylinder length: 


CL 



%clearance" 
“ ~ 100 


strk 
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Rate of molar storage for the control volume: 


1) 


dN 


cv 


dt 


with total number of moles, N 


l/out 


= JdN cv + N cv 

J ic 


cv 


Similarly, the rate of heavy gas molar storage is: 
dNhg cu V. v V , 


dt 


2 


F. X. 

in in 


F OUt X CV 


with total moles of heavy gas given as: Nhg cv = I dNhg cv +Nhg 


J. 


*CV. 
IC 


The heavy gas concentration in the control volume is then given as: 


Nhg 


X = 
cv N 


cv 


cv 


Rate of energy storage for the control volume: 
dE„ 


2) 


"CV 


dt 


Ye, F. T. - 'Sj 

X-U p in m m £-J 


C P cv F out T ou, + ^CV " W cv 


with, W rv = P 


dV 


CV 


CV CV dt 


The derivative of the control volume is evaluated as: 


3) 


dV 


CV 


^strk^ 


dt 


v 2 , 


w cos (wt + <>) A, 


cv 


It is assumed the gas behaves as an ideal gas with: 
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P„„ = N 


RT 


cv 


cv v 


cv 


Differentiating: 


dP 


dT 


dN 


N T dV 

'xv 'cv u V cv 


CV It N MI CV j "'’cv cv cv cv 

J dt " v^[ cv dt cv dt ' V cv dt 
Differentiate the energy equation, E C v = N cv Cv T cv , and rearrange to get: 


5) 


dT 


cv 


dt 


1 


dE 


cv 


c w dt 


- T, 


dN 


cv 


cv 


dt 


N 


cv 


Equations 1 , 2 and 3 along with the equation of state are substituted into 
Equations 4 and 5 which are than integrated to determine the system temper- 
ature and pressure with time. For fixed boundary volumes, dVcv/dt is set to 0. 


2.2 Branch Flow Analysis 


Flow through a restriction or frictional loss is represented by an incompressible 
resistance equation of the form: 

2 

AP = KpV , V = gas velocity 


Which can be written as: 


AP 



Wt 


F 2 


2.3 Evaluation of Heavy Gas Condensate Flow 

Molar air flow through the heat exchanger remains constant while heavy gas 
flow in the process stream is reduced by the amount of material condensed. For 
a constant temperature heat exchanger, condensate flow may be determined 
from inlet and outlet concentrations which are known. Figure 3 shows a 
schematic of the heat exchanger. 
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Heavy gas concentration is defined as: 


x 



( F . + F h ) 

^ air hgj 


with F . being constant for the process 


Entering molar flow of Heavy gas is given as: 



Similarly, exiting molar flow of heavy gas is given as: 

X„ 


hg2 


F . 

air ' 1 - X. 


Condensate flow is the difference in inlet and outlet heavy gas molar flow: 


F = F. - F = F . 

cond hgl hg2 air 


(i-x,)(i-x 2 ) J 


2.4 Valve Dynamics and Control 


Control valves in the system are pneumatically positioned. The positioner and 
valve dynamics together are taken as a first order system. The time constant for 
the valve dynamic system was either measured or assumed. 


Valve response is given by: Y do{ 



Tau 


This equation is integrated to give valve position: 


■ j Y do, dt + 


ic limit yL 


The control system is simple proportional integral (PI) control where the 
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error between the control variable and process variable is: err = P - P S et 


Proportional and integral commands are: 

P cmd = K p err 

f 

I = iK.errdt + ic limit I . L for anti-wind up integral control 
cmd J i cmd 0 r 

with command to the valve given as: Y C md = Pcmd + land 

2.5 Relief Valve Analysis 

The compressor system has relief valves located approximately at the heat 
exchanger for each stage. These relief valves protect the heat exchangers, 
which have the lowest rated design pressure in the system. The model 
evaluates system overpressure due to volume/relief valve interaction. 

For the relief valves, compressible nozzle flow analysis is used with mass flow 
given as: 


m = ^A* M a J yg c RTM wt 
a* 

For choked flow, M a = 1 , — = 1 

For an ideal gas, the ratio of static to stagnation temperature and pressure can 
be related by specific heat ratios: 

l 

, xY-1 

_2_J X = 2 

j+IJ ’ T t Y+1 

Substituting into the flow equation and solving for molar flow gives: 
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This equation agrees very well with the manufacturer's published data of flow 
versus pressure for these relief valves. 

3.0 Method of Solution 

A computer model of the above equations (Section 7.0) was generated using 
the ACSL code. ACSL is written to facilitate the solution of a system of nth 
order, nonlinear ordinary differential equations which describe a physical 
system. Special functions are available in addition to standard FORTRAN 
commands to provide flexibility in the problem description. The code sorts the 
commands into an executable sequence such that a variable is calculated 
before it is used, allowing the system of equations to be inputted in a logical 
sequence. Nine integration algorithms are available (default is fourth order 
Runga-Kutta) in addition to a user supplied routine. The code assigns lower 
order derivatives a variable name, then simultaneously solves the system of first 
order equations and determines the system states. The code allows input of 
user defined transfer functions for simplified plant or controller simulation. 

ACSL has extensive graphical capability and allows plotting of data from 
external files. 

4.0 Simulations 

Simulations of the compressor system were performed to validate the 
mathematical model and to evaluate compressor system operating conditions. 
Pressure transducers located in the suction and discharge side of each stage 
provide a pressure history of the compressor. Data was compiled for various 
compressor operations and compared with predicted data. 

The operator of the compressor system controls the suction pressure to the 
machine, discharge pressure of the machine, and 3rd stage discharge 
pressure. A concentration analyzer monitors the fraction of heavy gas at the 
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machine suction. The pressure and temperature data of the machine and 
heavy gas analyzer data is recorded on a PC at a rate of six data sets per 
minute. 

4.1 Steady State Analysis 

Simulated steady state conditions were initially compared with Dresser-Rand 
(Reference 2) calculated data and compressor measured data. Several cases 
of air mode calculated data from Reference 2 were compared with simulated 
data and were found to agree within 5%. Compressor data collected March 9, 
1995 was compared with calculated data. The data from this run includes a 
steady state condition with air mode, then a transition in inlet concentration and 
discharge pressure. The measured steady state air mode pressure data was 
compared with calculated data and is shown in Table 1 . Slight adjustments in 
the interstage pressure loss coefficient were made for better agreement. As can 
be seen from the table, there is good agreement between calculated and 
measured steady state data. In addition to the above condition, several steady 
state cases were run. Table 2 shows calculated conditions typical of nominal 
compressor operation including effects of inlet concentration variation. The 
results show good agreement with observed operation. Note that for high R12 
concentrations, condensation in the 4th stage goes to zero. This occurs due to 
the bypass flow from the 5th stage being at 21% FI 2 mixed with the 66% 
concentration of the 3rd stage discharge. These two gaseous streams mix for a 
concentration of 30% which is less than the R12 partial pressure at the 4th 
stage condensing temperature (43% partial pressure). Also note from the table 
that the inlet molar flow remains fairly constant with inlet concentration. 
However, the inlet mass flow increases with molecular weight which is reflected 
in higher interstage pressure drop. 

Figure 4 was constructed to show the effects of increased inlet pressure for a 
constant discharge pressure of 415 psia. Note that above about 6 psia suction 
pressure, the 4th stage discharge pressure exceeds the 5th stage discharge 
pressure by the amount of interstage pressure drop. 
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Table 1. Steady State Data Comparison: Air Mode Data of 3/9/95 


Location 

measured 

calculated 

difference 

measured 

calculated 

difference* 


(press psia) 

(press psia) 

(%) 

(temp °F) 

(temp °F) 

(%) 

I s * stage suction 

2.75 

2.75 

code input 

constant ~ 80 °F 


1 s * stage discharge 

10.12 

10.17 

0.49 

293.30 

316.04 

3.02 

2 nc * stage suction 

7.43 

7.59 

2.15 

constant ~ 80 °F 


2 nc * stage discharge 

36.39 

34.92 

4.04 

382.09 

362.84 

2.29 

3 f d stage suction 

32.60 

31.11 

4.57 

constant ~80 °F 


3 r d stage discharge 

204.47 

204.47 

code input 

462.06 

450.75 

1.23 

4^ stage suction 

190.86 

191.11 

0.13 

constant -80 °F 


4^ stage discharge 

380.26 

383.27 

0.79 

204.38 

191.54 

1.93 

5™ stage suciton 

368.49 

371.29 

0.76 

constant -80 °F 


5 th stage discharge 

465.61 

465.61 

code input 

126.47 

112.76 

2.34 


* Base on absolute (Rankine) scale 



Table 2. Calculated Steady State Conditions for the Clark Compressor 


inlet 

concentration 

0. 

21% 

33% 

50% 

66% 

95% 

Plsuc 

Pldis 

5.5 

22.07 

5.5 

23.33 

5.5 

23.93 

5.5 

24.59 

5.5 

25.13 

8.00 

34.24 

P2suc 

P2dis 

17.27 

69.76 

16.90 

72.13 

16.79 

73.89 

16.61 

75.17 

16.47 

76.30 

20.74 

97.25 

P3suc 

P3dis 

52.32 

200.14 

48.95 

194.39 

48.11 

200.00 

46.51 

200.00 

45.35 

200.00 

52.33 

200.00 

P4suc 

P4dis 

179.11 

396.2 

164.67 

385.67 

164.83 

390.43 

158.33 
351 .45 

152.61 

312.82 

143.90 

375.93 

P5suc 

P5dis 

375.47 

615.02 

355.46 

615.00 

354.42 

615.01 

314.96 

615.00 

278.89 

615.00 

320.30 

615.01 

suction 

flow 

0.129 

0.1138 

0.1075 

0.1006 

0.0953 

0.1119 

f"cond3 

0. 

0. 

0. 

0. 

0.0002 

0.1024 

f“cond4 

0. 

0. 

0. 

0.0177 

0.0343 

0. 

Fcond5 

0. 

0. 

0.0156 

0.0183 

0.0187 

0.0095 

bypas 

flow 

0. 

0. 

0.00402 

0.0088 

0.0126 

0.07732 

Xin 

0. 

0.21 

0.33 

0.50 

0.66 

0.95 

^3rd 

0. 

0.21 

0.33 

0.50 

0.659 

0.6593 

^4th inlet 

0. 

0.21 

0.325 

0.476 

0.606 

0.2964 

^4th 

0. 

0.21 

0.325 

0.375 

0.4215 

0.2964 

Xsth 

0. 

0.21 

0.214 

0.2144 

0.2144 

0.2144 

discharge flow 

0.127 

0.1136 

0.0915 

0.0636 

0.0406 

0.007857 


All pressures in psia, all flows in moles/sec. 

Fcond3,4,5 - flows from collection pot. 

bypas flow • flow from 5th stage discharge to 4th stage suction 

When comparing with measured data, suction pressures should be the same for valid comparison. 
Mwt air = 28.97 
Mwt FI 2 = 120.93 

Mass flow of collection pots = (moles/sec)*MwtFl2 = Ibm/sec 

Mass flow suction, discharge, and bypass = (moles/sec)*(X*M w tFl 2 +(‘l'X)*M w tair) = Ibm/sec 
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4.2 Transient Analysis 


A comparison was made of measured data versus calculated for an operational 
run conducted March 9, 1995. For this run, approximate steady conditions were 
achieved with very low inlet concentration of about 1% R12, suction pressure of 
2.74 psia, 3rd stage discharge of 205 psia, and a machine discharge pressure 
of 470 psia. From the steady conditions, inlet pressure was reduced, then 
increased, held steady, and then increased again. During this period, inlet 
concentration varied from about 1% up to 92% R12. The 3rd stage discharge 
pressure and machine discharge pressure were also varied during this period. 
This transition period covered about 4 V 2 minutes during the run of March 9, 
1995. It should be noted that the suction pressure, 3rd stage discharge 
pressure, and machine back pressure set points are controlled manually by the 
operator from the control panel. These set points are dialed in by the operator 
to achieve conditions deemed optimum. The set points are not recorded; 
however, the resultant process variable is. In order to assess the math model 
transient performance, the input process variables of inlet concentration, suction 
pressure, 3rd stage discharge pressure, and machine back pressure were used 
as set points to the math model. These set points were compared to process 
variables and the error used to generate a command to the control valve as 
discussed in Section 2. Figures 5 through 15 compare calculated versus 
measured data for inlet concentration, stage suction and stage discharge 
pressures. Recorded data is shown as triangles and calculated data as solid 
lines. 

Comparison of the calculated data with measured data shows agreement within 
about 20% for variable magnitudes. Some error between the calculated and 
measured data is due to the method of introducing control set points. Figures 5, 
6, 11, and 15 show the control variables. The math model process variables 
are a result of using measured process variables as set points as opposed to 
the actual set points; therefore, the differences in control gain modeling can 
manifest in differences in process variables between measured and calculated. 
It is, however, concluded that the math model transient response is 
representative of the actual machine and, therefore, is suitable for use in the 
prediction of machine response to transient conditions. 
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4.3 Relief Valve Dynamic Analysis 


An analysis of the relief valves and their dynamic interaction with the 
compressor system was performed for several cases. Case 1 evaluated the 
relief valve capacity relative to worse case pressure condition resulting from a 
failure of the discharge valve to the closed position. From steady state 
conditions of 5.5 psia suction and 615 psia discharge, the back pressure valve 
was commanded shut at 80 seconds into the simulation. Figure 16 shows valve 
position with time. The valve completely shuts in about 4 seconds. The 
discharge pressure, shown in Figure 17, increases as flow continues to fill the 
discharge volume. Approximately 7 seconds into the event, the 5th stage 
pressure reached the relief valve cracking pressure of 720 psia and flow 
through the relief valve was established as shown in Figure 18. The results of 
this simulation show that the relief valve capacity is adequate for this type of 
failure. 

As stated previously, R134a refrigerant is combustible under certain conditions. 
From the combustibility tests, a molar heat of combustion for a stoichiometric 
mixture and approximate flame velocity were determined for a gas/air mixture. 1 
The burning dynamics were modeled with ignition occurring at the cylinder 
discharge. Five simulations were evaluated for a worst case ignition scenario 
for each stage. From a steady state condition of 5.5 psia suction and 615 psia 
discharge, the stage flow entering the adjacent upstream stage, or discharge 
valve for the 5th stage, was set to zero. Simultaneously ignition was initiated at 
the cylinder discharge. Pressure rises in the trapped gas as combustion 
proceeds and exit flow is held to zero. Results for the 4th stage analysis are 
shown in Figures 19 and 20. From steady state conditions, ignition was initiated 
at 50 seconds, simultaneously setting volume discharge flow to zero. Figure 19 
shows 4th stage discharge pressure increase due to combustion and blocked 
flow to reach approximate steady conditions. The relief valve cracks about 3 V 2 
seconds into the event. It's molar flow is shown in Figure 20. Results for all 
cases are shown in Table 3. 


1 Babcock, D. A., Bruce, R. A., "Combustibility Tests of 1,1,1,2-tetrafluoroethane in a 
Simulated Compressor Cylinder." Unpublished Report, 1995. 
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Table 3. Relief Valve Capacity and Stage Overpressure Due To Combustion 


stage 

1st 

2nd 

crack press 
(psig) 

(psia) 

35 

49.7 

105 

119.7 

capacity* 

(SCFM) 

25686 

24643 

3537 

peak pressure 
due to comb, (psia) 

52 

169 

max allowable 
pressure** (psia) 

290 

620 

peak relief valve 
flow (moles/sec) 
due to comb. 

0.175 

0.29 


3rd 

4th 

5th 

350 

600 

705 

364.7 

614.7 

719.7 

13089 

22094 



450 

652 

780 

1620 

1960 

1960 

0.30 

0.50 

0.55 


Capacity at 60 deg F and 10% overpressure 
^internal pressure at ultimate 



4.4 Control System Evaluation 


The control system consists of operator set PI pressure control for the suction 
pressure, the third stage discharge pressure, and machine discharge pressure. 
Consideration is being given to scheduling third stage pressure and discharge 
pressure set points as a function of machine suction pressure and inlet 
concentration. A scheme is implemented to maintain a pressure ratio of 6 or 
less across any individual stage of the compressor system. That scheme will be 
considered here. 

As discussed, the lower limit for auto-ignition of a R134a/air mixture is about 
500 °F. 1 Such a condition may occur if the isentropic compression ratio for any 
cylinder exceeds 6.6 with an inlet temperature of 100 °F and specific heat ratio 
of 1 .4 (air). It should be noted that as the concentration of heavy gas increases, 
specific heat ratio decreases, greatly increasing required pressure ratio for a 
500 °F discharge temperature. It was decided, however, to use a pressure ratio 
set point of 6.0 for each cylinder regardless of concentration. Referring to 
Figure 1 , the pressure ratio for each cylinder is calculated and compared to the 
set point. If the cylinder pressure ratio is less than the set point, the bypass 
valve from 5th stage discharge to machine suction remains closed. If the 
cylinder pressure ratio exceeds the set point, the PI controller generates a 
command to the valve which diverts flow from the machine discharge to suction, 
proportional to the error. This method of control raises the suction pressure and 
lowers the discharge pressure such that machine and individual stage pressure 
ratios are reduced only by that required to maintain a ratio of 6 or less. As the 
transient or abnormal operating condition is alleviated, normal control and 
operation is restablished. 

An example using the model to predict compressor performance is presented 
for a scenario in which the suction valve is commanded shut from a normal 
steady state operating condition. With the simulation at steady state conditions 
of 5.5 psia suction and 615 psia discharge, the suction pressure set point is set 
to zero 65 seconds into the simulation. Such an action could occur due to 
operator error or failure of the suction control valve. Figure 21 shows the action 

1 Babcock, D. A., Bruce, R. A., "Combustibility Tests of 1,1,1,2-tetrafluoroethane in a 
Simulated Compressor Cylinder." Unpublished Report, 1995. 
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of the suction control valve. Figure 22 shows the decrease with time of the 
suction pressure. Due to the decrease in suction pressure, the pressure ratio 
(pressure ratio = discharge/suction) across each stage increases with time. The 
pressure ratio across the 2nd stage is shown in Figure 23. The 2nd stage 
pressure ratio reaches a value of 6 about 1 5 seconds into the event and 
reaches this value before the other stages. 

This same simulation was then run with a pressure ratio set point of 6 for any 
stage. The second stage reaches the limiting ratio first and regulation begins. 
Figure 24 shows the 2nd stage pressure ratio versus time. As before, the upset 
begins 65 seconds into the simulation. At about 80 seconds into the simulation, 
the 2nd stage pressure ratio reaches 6 and regulation begins. From Figure 24, 
it is seen there is some overshoot and then oscillation in stage pressure ratio as 
regulation begins. Molar flow from the discharge to the suction side of the 
machine is shown in Figure 25. The flow appears to oscillate at about .02 
moles per second. This flow serves to raise the suction pressure and, therefore, 
reduce the pressure ratio across each stage. Figure 26 shows machine suction 
pressure for the controlled case and is compared to Figure 22 which is the 
uncontrolled suction pressure. 

Results of this simulation show that such a control scheme is viable. The model 
allows controller gain evaluation and simulation of various failure or abnormal 
operating scenarios. 

5.0 Conclusion 

A mathematical model of a multi-stage compressor with variable molecular 
weight flow medium was derived. The equations of state were solved using the 
Advanced Continuous Simulation Language (ACSL). Results of simulated 
operation were compared with measured data. The code has input variables 
and constants defined initially and can be changed interactively to easily 
assess the impact of hardware changes or changes in control scheme. The 
code is in a general format making application to other multi-stage compressor 
systems simple. 
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8.0 Code Listing 

program dark compressor 

"234567890123456789212345678931234567894123456789512345678961234567897 


TDT CLARK COMPRESSOR PERFORMANCE 
ANALYSIS 

D. A. BABCOCK JUNE 1995 


nomenclature 


f 

fdt 

flow 

rate of change 

of 

flow 

moles/sec 

moles/sec**^2 

n 

ndt 

total number of moles 
rate of change of moles 

n 

n/sec 

e 

edt 

energy 

rate of change 

of 

energy 

btu 

btu/sec 

tc 

tcdt 

temperature 
rate of change 

of 

temperature 

deg R 
deg R/sec 

P 

pdt 

pressure 
rate of change 

of 

pressure 

psi 

psi/sec 


constant cpa 
constant cva 
constant r 
constant ma 

it 

M 

It 

It 


It 

tt 

tl 

constant alsl 
constant als2 
constant als3 
constant als4 
constant a2sl 
constant a2s2 
constant a3sl 
constant a3s2 
constant a4sl 
constant a4s2 
constant a5sl 
constant a5s2 


= 7,1846 

= 5.1277 

= 1545.0 

= 28.97 


24 . 68 

22.69 
120.93 

131.86 

3761.0 


= 7.0686 

= 7.0195 

= 7.0686 

= 7.0195 

= 4.2761 

= 4.2270 

= 1.3104 

= 1.2613 

= 0.3712 

= 0.3221 

= 0.1963 

= 0*. 1473 


constant cphg = 
constant cvhg = 
constant mhg = 
constant vplOO = 
constant hliq = 

it 


constants for air 

$"btu/ (lbmole-deg r) , const pres spec ht 
$"btu/ ( lbmole-deg r) # const vol spec ht 
$”ft-lbf/ (lbmole-deg r) , ideal gas const 
$"lbm/lbmole, molecular wt air 


constants for R-12 

$"btu/ (lbmole-deg r) , const pres spec ht 
$"btu/ (lbmole-deg r) , const vol spec ht 
$"lbm/lbmole, molecular wt R-12 
$"psia, sat vapor at 100 deg F 
$"btu/lbmole, liquid enthalpy 


geometry constants 

$"ft**2, cross section area stage 1 cyl 
$"stage 1 cross section less rod area 
$"f t**2 , cross section area stage 1 cyl 
$"stage 1 cross section less rod area 
$"f t**2 , cross section area stage 2 cyl 
$"stage 2 cross section less rod area 
$ H ft**2, cross section area stage 3 cyl 
$"stage 3 cross section less rod area 
$”f t**2 , cross section area stage 4 cyl 
$ M stage 4 cross section less rod area 
$”f t**2 , cross section area stage 5 cyl 
$"stage 5 cross section less rod area 



constant 

strk = 

1.4167 

constant 

vefl = 

1.13 

constant 

vef 2 = 

1.1 

constant 

vef 3 = 

1.06 

constant 

vef 4 = 

1.03 

constant 

•i 

vef 5 = 

1.02 

constant 

V0 = 

3284.0 

constant 

vl = 

315.0 

constant 

v2 = 

30$. 0 

constant 

v3 = 

208.0 

constant 

v4 

146.0 

constant 

V5 = 

63.0 

constant 

H 

v3m = 

10.0 

If 

|| 


other 

constant 

pcbl=4 15 . 0 

constant 

cg=66 . 

0 

constant 

tstrt 

= 0.0, 

M 


control 

constant 

set 

= 615.0 

constant 

setO 

= 5.5 

constant 

prset 

= 20.0 

constant 

setb 

= 0.0 


$"ft, piston stroke 
^"compression eff 
$"compression eff 
$"compression eff 
$" compress ion eff 
$"compression eff 

$"f t**3 , suction volume 
$ " f t * * 3 , dwn stream first stage 
$"ft**3 , dwn stream second stage 
$"f t**3 , dwn stream third stage 
$»»ft**3, dwn stream fourth stage 
$"ft**3, dwn stream fifth stage 
$"f t**3 , dwn stream fifth stage 


constant 

constant 

constant 

m 

constant 

constant 

constant 

constant 

tt 

constant 

constant 

constant 

constant 

it 

constant 

constant 

constant 

constant 


$"psia, pressure dwnstrm of back pres vlv " 

$"back pressure valve flow coeficient 11 

tstp = 1.0 

ler set points ' 

$"psia, pressure set point machine disch * 

$"psia, pressure set point suction press * 

$"psia, pressure set point pratio set point' 
$"psia, pressure set point 3rd stg dis 
discharge controller 1 

kp = 0.05 $"proportional constant 1 

ki = .01 $"integral constant • 

tau = 0.95 $"sec, valve time constant 1 

ratio controller 1 

kpl = 0.1 $"proportional constant 1 

kil = .05 $"integral constant 1 

taul = 1.0 $"sec, inlet valve time constant 1 

vcnstl = 12.0 $"psi, constant for flow cntrl valve 1 

suction controller 1 

kpO =0.13 
kiO = .055 
tauO = 0.75 
vcnst = 800.0 


$ "proportional constant 
$"integral constant 
$"sec, inlet valve time constant 
$"psi, constant for flow cntrl valve 


3rd stage pressure controller 


kpb = 0.05 
kib = .001 
taub = 1.00 
vcnstb= 800.0 


$"proportional constant 
$" integral constant 
$"sec, inlet valve time constant 
$"psi, constant for flow cntrl valve 


other constants 


constant j = 778.16 $"f t-lbf /btu, conversion factor 

constant omega = 31.415927 $"rad/sec, 300 rpm machine 
constant pi = 3.1415927 $"rad, constant 

i« 

" relief valve constants 


constant al = 25.25 , a2 = 10 . 734 , a3 = l . 784 , a4=l - 784 , a5=l . 784 
" these are the choked areas for the relief valves 


constant cpl=49.7 
constant cp2=164.7 
constant cp3=364.7 
constant cp4=614.7 


, fpl=53 . 2 
, fp2=179 . 7 
, f p3=399 . 7 
, f p4 = 674 • 7 
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constant cp5=719.7 , fp5=790.2 

constant cnst=63.69 

" cp is the crack pressure for the relief valve, fp is the full flow " 
" pressure of the relief valve " 

•i •' 

" following are the loss coefficients for the compressor valves " 

constant kvO=. 15588 

constant klslup=l . 5 , kls2up=l . 5 , kls3up=l . 5 , kls4up=l . 5 

constant klsldn=l . 5 , kls2dn=l . 5 , kls3dn=l . 5 , kls4dn=l . 5 

constant k2slup=2 . 03 , k2s2up=2 . 03 

constant k2sldn=2.03,k2s2dn=2.03 

constant k3slup=2 . 9 ,k3s2up=2 . 9 

constant k3sldn=2.9,k3s2dn=2.9 

constant k4slup=2.9,k4s2up=2.9 

constant k4sldn=2.9,k4s2dn=2.9 

constant k5slup=2.9,k5s2up=2.9 

constant k5sldn=2 . 9 ,k5s2dn=2 . 9 


initial $ M initial conditions of state variables 1 

n ' 

losl2 = 0.23 
los23 = 2.86 
los34 = 2.1014 
los45 = 19.71 

kv3m = 7.434 $"loss coef from vol 3 to vol3 mix chamber 

ii 

» initial temperatures and pressures of fixed volumes 

II 

tv3mic = 560.0 $ tv3m = tv3mic $ pv3mic = 14.7 $ pv3m=pv3mic 
tsuc = 560.0 $ psuc =14.7 $ ycmd0=.5 $ y0=.5 

tvOic = 560.0 $ tvO = tvOic $ pvOic =14.7 $ pvO = pvOic 

tvlic = 560.0 $ tvl = tvlic $ pvlic =14.7 $ pvl = pvlic 

tv2ic = 560.0 $ tv2 = tv2ic $ pv2ic = 14.7 $ pv2 = pv2ic 

tv3ic = 560.0 $ tv3 = tv3ic $ pv3ic =14.7 $ pv3 = pv3ic 

tv4ic = 560.0 $ tv4 = tv4ic $ pv4ic =14.7 $ pv4 = pv4ic 

tv5ic = 560.0 $ tv5 = tv5ic $ pv5ic =14.7 $ pv5 = pv5ic 

ii » 

pf slic=pv0ic $ pfs2ic=pvlic $ pfs3ic=pv2ic $ pfs4ic=pv3ic 
pf s5ic=pv4ic $ pfdlic=pvlic $ pfd2ic=pv2ic $ pfd3ic=pv3ic 
pfd4ic=pv4ic $ pfd5ic=pv5ic 

tfdlic=tvlic $ tfd2ic=tv2ic $ tfd3ic=tv3ic $ tfd4ic=tv4ic 
tfd5ic=tv5ic 

" 

" initial temperatures and pressures of cylinder volumes " 

" dark 1st stage " 

tslic=100. $ ts2ic=100. $ ts3ic=100. $ ts4ic=100. $ ts5ic=100. 
tsl=tslic $ ts2=ts2ic $ ts3=ts3ic $ ts4=ts4ic $ ts5=ts5ic 
tl=tsl $ t2=ts2 $ t3=ts3 $ t4=ts4 $ t5=ts5 

tlslic = 560.0 $ tlsl = tlslic $ plslic = 14.7 $ plsl = plslic 

tls2ic = 560.0 $ tls2 = tls2ic $ pls2ic = 14.7 $ pls2 = pls2ic 

tls3ic = 560.0 $ tls3 = tls3ic $ pls3ic = 14.7 $ pls3 = pls3ic 

tls4ic = 560.0 $ tls4 = tls4ic $ pls4ic = 14.7 $ pls4 = pls4ic 

" dark 2nd stage " 

t2slic = 560.0 $ t2sl = t2slic $ p2slic = 14.7 $ p2sl = p2slic 

t2s2ic = 560.0 $ t2s2 = t2s2ic $ p2s2ic = 14.7 $ p2s2 = p2s2ic 

" dark 3rd stage 

t3slic = 560.0 $ t3sl = t3slic $ p3slic = 14.7 $ p3sl = p3slic 

t3s2ic = 560.0 $ t3s2 = t3s2ic $ p3s2ic = 14.7 $ p3s2 = p3s2ic 

" dark 4th stage 

t4slic = 560.0 $ t4sl = t4slic $ p4slic = 14.7 $ p4sl = p4slic 
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t4s2ic = 560.0 $ t4s2 = t4s2ic $ p4s2ic = 14.7 $ p4s2 = p4s2ic 
" dark 5th stage 

tSslic = 560.0 $ t5sl = t5slic $ p5slic = 14.7 $ p5sl = p5slic 
t5s2ic = 560.0 $ t5s2 = t5s2ic $ p5s2ic = 14.7 $ p5s2 = p5s2ic 

•i 

tin = 540.0 $ M all cylinder inlet temperatures are 100 deg f 
xin =0.0 Centering concentration of HG, hg/air moler basis 

t» initial moles of HG 

nhgOic = pv0ic*144 , *v0/ (r*tv0ic) *x0 
nhglic = pvlic*144 . *vl/ (r*tvlic) *xl 
nhg2ic = pv2ic*144 . *v2/ (r*tv2ic) *x2 
nhg3ic = pv3ic*144 . *v3/ (r*tv3ic) *x3v 
nhg4ic = pv4ic*144 . *v4/ (r*tv4ic) *x4v 
nhgSic = pv4ic*144 . *v5/ (r*tv5ic) *x5v 

it 

" compressor efficiency 

clsl=strk* ( . 5+ . 078*vef 1) 

cls2=strk* ( . 5+ . 076*vef 2 ) 

cls3=strk* ( . 54*. 070*vef 3 ) 

cls4=strk* ( . 54-. 126*vef4) 

cls5=strk* ( . 54- . 3 4 0*vef 5 ) 

end $"of initial section 

dynamic 

derivative 


” volume 0 analysis 

tt 

procedural 

cpO = x0*cphg 4- (l.-x0)*cpa 
cvO = xO*cvhg 4- (l.-xO)*cva 
mO = xO*mhg 4- (l.-xO)*ma 
cpin = xin*cphg 4- (l.-xin)*cpa 
end 

it 

fO = sqrt (abs (psuc-pvO) ) * (psuc/tsuc) /kvO*sign (1 . 0 , psuc-pvO) *y0 

ydotO=(ycmdO-yO) * tauO 

yO = limint (ydotO, 0. 0, 0. 0 , 1 . 0) 

errO = -(pvO-setO) 

pcmdO=kpO*errO 

icmdO=limint (kiO*errO ,0.0,0. 0,1.0) 
ycmdO=bound ( 0 . 0 , 1 . 0 , pcmd04-icmd0 ) 

nvOdt = fO - (f lslup+f Is2up4-f ls3up+f ls4up) 4-f5to0 
nvO = pv0*144 . 0*v0/ (r*tvO) 

nhgOdt=f 0*xin4-f 5to0*x5- (f lslup4-f ls2up+f Is3up4-f ls4up) *x0 

nhgO = integ (nhgOdt, nhgOic 

x0=nhg0.nv0 

tvOdt = (evOdt/cvO-nvOdt*tvO) /nvO 
tvO = integ(tvOdt , tvOic) 

pvOdt = r/ (v0*144 . ) * (nvO*tvOdt+tvO*nvOdt) 
pvO = integ (pvOdt ,pvOic) 

evOdt = cpO* (fO*tsuc-(f lslup4*f Is2up4-f ls3up+f ls4up) *tvO) 4-cp5*f 5to0*tin 

it 

ii dark 1st stage analysis 

it 

flslup=rsw(psucl-plsl. le. 0. ,0. , sqrt (psucl-plsl) /klslup) 
f Is2up=rsw(psucl-pls2. le. 0. , 0. , sqrt (psucl-pls2) /kls2up) 
fls3up=rsw(psucl-pls3 . le. 0. , 0. , sqrt (psucl-pls3) /kls3up) 
fls4up=rsw(psucl-pls4 . le. 0. , 0. , sqrt (psucl-pls4) /kls4up) 
nlsldt = f lslup-f lsldn 
nls2dt = f ls2up-f ls2dn 



nls3dt = f ls3up-f ls3dn 

nls4dt = f ls4up-f ls4dn 

nisi = plsl*144.*vlsl/(r*tlsl) 

nls2 = pls2*144 . *vls2/ (r*tls2) 

nls3 = pls3*144.*vls3/(r*tls3) 

nls4 = pls4*144 . *vls4/ (r*tls4) 

tlsldt=(elsldt/cvO-nlsldt*tlsl) /nisi 

tls2dt= (els2dt/cvO-nls2dt*tls2 ) /nls2 

tls3dt=(els3dt/cvO-rnls3dt*tls3) /nls3 

tls4dt=(els4dt/cvO-nls4dt*tls4) /nls4 

tlsl = integ(tlsldt, tlslic) 

tls2 = integ(tls2dt , tls2ic) 

tls3 = integ(tls3dt, tls3ic) 

tls4 = integ(tls4dt, tls4ic) 

plsldt=r/vlsl* (nlsl*tlsldt+tlsl*nlsldt-nlsl*tlsl/vlsl*vlsldt) /144 . 
pls2dt=r/vls2* (nls2*tls2dt+tls2*nls2dt-nls2*tls2/vls2*vls2dt) /144 . 
pls3dt=r/vls3* (nls3*tls3dt+tls3*nls3dt-nls3*tls3 /vls3 *vls3dt ) / 14 4 . 
pls4dt=r/vls4* (nls4*tls4dt+tls4*nls4dt-nls4*tls4/vls4*vls4dt) /144 . 
plsl = integ(plsldt,plslic) 
pls2 = integ(pls2dt ,pls2ic) 
pls3 = integ(pls3dt,pls3ic) 
pls4 = integ(pls4dt ,pls4ic) 

it 

vlsl=(clsl*vef l-strk/2 . *sin (omega*t) ) *alsl 
vls2=(clsl*vef l-strk/2 . *sin(omega*t+pi) ) *als2 
vis 3= ( els l*vef l-strk/2 . *sin (omega*t+pi) ) *als3 
vls4= (clsl*vef l-strk/2 . *s in (omega *t) ) *als4 

it 

vlsldt=-alsl*strk/2 . *omega*cos (omega*t) 
vls2dt=-als2*strk/2 . *omega*cos (omega*t+pi) 
vls3dt=-als3*strk/2 . *omega*cos (omega*t+pi) 
vls4dt=-als4*strk/2 . *omega*cos (omega*t) 

ti 

elsldt=cpO* (f lslup*tin-f lsldn*tlsl) -plsl*vlsldt*144 . / j 
els2dt=cp0* (f ls2up*tin-f Is2dn*tls2) -pls2*vls2dt*144 . / j 
els3dt=cp0* (f ls3up*tin-f Is3dn*tls3) -pls3*vls3dt*144 . / j 
els4dt=cpO* (f ls4up*tin-f Is4dn*tls4) -pls4*vls4dt*144 . / j 

II 

tl=rsw(f lsldn. gt . 0. ,tlsl-460. , tl) 

ti 

" volume 1 analysis 

it 

f lsldn=rsw(plsl-pvl . le. 0. , 0. , sqrt (plsl-pvl) /klsldn) 
fls2dn=rsw(pls2-pvl. le.O. ,0. , sqrt (pls2-pvl) /kls2dn) 
f Is3dn=rsw(pls3-pvl. le.O. ,0. , sqrt (pls3-pvl) /kls3dn) 
f Is4dn=rsw(pls4-pvl. le.O. ,0. , sqrt (pls4-pvl) /kls4dn) 

II 

nvldt = (f lsldn+fls2dn+f ls3dn+f ls4dn) - ( f 2slup+f 2s2up) 

nvl = pvl*144.*vl/ (r*tvl) 

nhgldt = flin*xO - flOut*xl 

nhgl = integ (nhgldt, nhglic) 

xl = nhgl/nvl 

tvldt = (evldt/cvO-nvldt*tvl) /nvl 
tvl = integ(tvldt, tvlic) 

pvldt = r/ (vl*144 . ) * (nvl*tvldt+tvl*nvldt) 
pvl = integ (pvldt, pvlic) 

it 

evldt=cp* (f lsldn*tlsl+f Is2dn*tls2+f Is3dn*tls3+f Is4dn*tls4) tenrout 
enrout = - (f2slup+f 2s2up+f lowl) *tvl*cp + qburnl 



tgol = rsw(tl.gt. 500.0. or. tgol.eq.l. 0,1. 0,0.0) 
qburnl=hcomb* (f Is ldn+f ls2dn+f ls3dn+f ls4dn+. 00054 *pvl) *tgol 
ratiol=(pvl-cpl) / (fpl-cpl) 

fulll=rsw(pvl-cpl. le. 0. ,0. , fpl*al/sqrt (tvl) /cnst) 
f lowl=rsw(ratiol. It. 1. O,ratiol*fulll,pvl/fpl*fulll) 

i» « 

H ” 

» dark 2nd stage analysis " 

ii « 

f2slup=rsw(psuc2-p2sl.le.0. ,0. f sqrt (psuc2-p2sl) /k2slup) 

f 2s2up=rsw(psuc2-p2s2 . le, 0 . , 0. , sqrt (psuc2-p2s2) /k2s2up) 

n2sldt = f 2slup-f2sldn 

n2s2dt = f 2s2up-f 2s2dn 

n2sl = p2sl*144 . *v2sl/ (r*t2sl) 

n2s2 = p2s2*144 . *v2s2/ (r*t2s2) 

t2sldt=(e2sldt/cv0-n2sldt*t2sl) /n2sl 

t2s2dt= (e2s2dt/cv0-n2s2dt*t2s2 ) /n2s2 

t2sl = integ(t2sldt, t2slic) 

t2s2 = integ (t2s2dt , t2s2ic) 

p2sldt=r/v2sl* (n2sl*t2sldt+t2sl*n2sldt-n2sl*t2sl/v2sl*v2sldt) /144 . 
p2s2dt=r/v2s2* (n2s2*t2s2dt+t2s2*n2s2dt-n2s2*t2s2 /v2s2*v2s2dt) /144 . 
p2sl = integ(p2sldt,p2slic) 
p2s2 = integ(p2s2dt,p2s2ic) 

ti 

v2sl=(cls2*vef2-strk/2 . *sin (omega* t) ) *a2sl 
v2s2=(cls2*vef 2-strk/2 . *sin (omega*t+pi) ) *a2s2 

ii 

v2sldt=-a2sl*strk/2 . *omega*cos (omega*t) 
v2s2dt=-a2s2*strk/2 . *omega*cos (omega*t+pi) 

ii 

e2sldt=cp0* (f 2slup*tin-f 2sldn*t2sl) -p2sl*v2sldt*144 . / j 
e2s2dt=cp0* (f 2s2up*tin-f 2s2dn*t2s2) -p2s2*v2s2dt*144 . / j 


t2=rsw(f2sldn.gt. 0, ,t2sl-460. ,t2) 

it 

it volume 2 analysis 

it 

f 2sldn=rsw(p2sl-pv2 . le. 0, , 0. , sqrt (p2sl-pv2) /k2sldn) 
f 2s2dn=rsw (p2s2-pv2 . le. 0 • , 0 . , sqrt (p2s2-pv2 ) /k2s2dn) 
m 


nv2dt - ( f 2sldn+f 2s2dn) - ( f 3slup+f 3s2up) 

nv2 = pv2*144 . *v2/ (r*tv2) 

nhg2dt = f2in*x0 - f20ut*x2 

nhg2 = integ(nhg2dt,nhg2ic) 

x2 = nhg2/nv2 

tv2dt = (ev2dt/cv0-nv2dt*tv2) /nv2 
tv2 = integ(tv2dt, tv2ic) 

pv2dt = r/ (v2*144 . ) * (nv2*tv2dt+tv2*nv2dt) 
pv2 = integ (pv2dt,pv2ic) 

ii 

ev2dt=cp* (f 2sldn*t2sl+f 2s2dn*t2s2- (f 3slup+f3s2up+f low2) *tv2) +qburn2 

If 


tgo2 = rsw(t2 .gt . 500 . 0. or . tgo2 . eq. 1 . 0, 1. 0, 0. 0) 
qburn2=hcomb* (f 2sldn+f 2s2dn+. 00032*pv2) *tgo2 
ratio2=(pv2-cp2) / (fp2-cp2) 

full2=rsw (pv2-cp2 . le. 0. , 0. , fp2*a2/sqrt (tv2) /cnst) 
f low2=rsw(ratio2 . It. 1. 0, ratio2*full2 ,pv2/fp2*full2) 


if 


ii 


dark 3rd stage analysis 


if 



f 3slup=rsw(psuc3-p3sl . le. 0. , 0. , sqrt (psuc3-p3sl) /k3slup) 

f3s2up=rsw(psuc3-p3s2.1e.0. ,0. , sqrt (psuc3-p3s2 ) /k3s2up) 

n3sldt = f 3slup-f 3sldn 

n3s2dt = f 3s2up-f 3s2dn 

n3sl = p3sl*144 . *v3sl/ (r*t3sl) 

n3s2 = p3s2*144 . *v3s2/ (r*t3s2) 

t3sldt=(e3sldt/cv0-n3sldt*t3sl) /n3sl 

t3s2dt=(e3s2dt/cv0-n3s2dt*t3s2) /n3s2 

t3sl = integ(t3sldt, t3slic) 

t3s2 = integ(t3s2dt‘, t3s2ic) 

p3sldt-r/v3sl* (n3sl*t3sldt+t3sl*n3sldt-n3sl*t3sl/v3sl*v3sldt) /144 . 
p3s2dt=r/v3s2* (n3s2*t3s2dt+t3s2*n3s2dt-n3s2*t3s2/v3s2*v3s2dt) /144 . 
p3sl = integ(p3sldt,p3slic) 
p3s2 = integ(p3s2dt,p3s2ic) 

« n 

v3sl=(cls3*vef 3-strk/2 . *sin (omega *t+pi) ) *a3sl 
v3s2=(cls3*vef 3-strk/2 . *sin (omega* t) ) *a3s2 

ii M 

v3sldt=-a3sl*strk/2 . *omega*cos (omega*t+pi) 
v3s2dt=-a3s2*strk/2 . *omega*cos (omega*t) 

it it 

e3sldt=cp0* ( f 3slup*tin-f 3sldn*t3sl) -p3sl*v3sldt*144 . / j 
e3s2dt=cp0* (f 3s2up*tin-f3s2dn*t3s2) -p3s2*v3s2dt*144 . / j 

•• it 

t3=rsw(f 3sldn.gt . 0. ,t3sl-460. f t3) 

h ii 

i» volume 3 analysis " 

H n 

f 3sldn=rsw(p3sl-pv3 . le. 0. , 0. , sqrt (p3sl-pv3) /k3sldn) 
f 3s2dn=rsw(p3s2-pv3 . le. 0 . , 0. , sqrt (p3s2-pv3 ) /k3s2dn) 

i» it 

nv3dt = (f 3sldn+f 3s2dn) - (f 3out) 
nv3 = pv3 *144 . *v3 / (r*tv3 ) 
nhg2dt = f2in*x0 - f20ut*x2 
nhg2 = integ (nhg2dt , nhg2ic) 
x2 = nhg2/nv2 

tv3dt = (ev3dt/cv0-nv3dt*tv3) /nv3 
tv3 = integ (tv3dt f tv3ic) 

pv3dt = r/ (v3*144 . ) * (nv3*tv3dt+tv3*nv3dt) 
pv3 = integ (pv3dt,pv3ic) 

n ii 

ev3dt=cp0* ( (f 3sldn*t3sl+f 3s2dn*t3s2) - (f3out) *tv3) +qburn3 

ii it 

tgo3 = rsw(t3 .gt . 500 . 0 . or . tgo3 . eq. 1 . 0, 1 . 0, 0. 0) 
qburn3=hcomb* ( f 3sldn+f 3s2dn+ . 000116*pv3 ) *tgo3 
ratio3=(pv3-cp3) / (fp3-cp3) 

full3=rsw(pv3-cp3 . le. 0. ,0. , fp3*a3/sqrt (tv3) /cnst) 
f Iow3=rsw(ratio3 . It. 1 . 0 , ratio3*full3 , pv3/fp3*full3 ) 


ii ii 

ii ii 

u conditions leaving 3rd stage after cooler " 

ii ii 

x3calc=vplOO/pv3 

x3=rsw(pv3 .gt. vplOO . and.x3calc. It . xO , x3calc, xO) 
dp3m=pv3-pv3m 


f 3out=realpl (0. 5 , sqrt (abs (dp3m) *pv3/ (tv3*kv3m*m0) ) *sign ( 1 . , dp3m) ,0.0) 
f3air=f3out* (l.-xO) 

f cond3=f 3air* (xO* ( 1 . -x3 ) -x3 * ( 1 . -xO) ) / ( 1 . -xO) / ( 1 . -x3 ) 
f3mix= f 3out - fcond3 

»• i» 
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11 


5th stage discharge to 4th stage suction valve 


it 


ii " 

M valve dynamics and controller, assume 1st order " 

ii •• 

ydotb = (ycmdb-yb) *taub 

yb=bound(0. ,1. , limint (ydotb, 0 . 0 , 0 . 0 , 1 . 0) ) 
fbypas = yb* (pv5-pv3m) /vcnstb 
errb=rsw (pv5 . It . pv3m , -pv3 , setb-pv3 ) 
pcmdb=kpb*errb 

icmdb=limint (kib*errb, 0 . 0 , 0 . 0 , 1 . 0 ) 
ycmdb=bound ( 0 . 0 , 1 . 0 , pcmdb+icmdb) 


ti it 

” volume 3 mixing chamber " 

ii it 


procedural 

x3m=rsw(f 3mix+fbypas. le. 0. , 0, , (f 3mix*x3+fbypas*x5) / (f 3mix+f bypas) ) 

cp3m=x3m*cphg + ( 1 . -x3m) *cpa 

cv3m=x3m*cvhg + (l.-x3m)*cva 

m3m=x3m*mhg + (l.-x3m)*ma 

end 

n n 

nv3mdt = f3mix - ( f 4slup+f 4s2up) 4- fbypas 
nv3m = pv3m*144 . *v3m/ (r*tv3m) 
tv3mdt=(ev3mdt/cv3m-nv3mdt*tv3m) /nv3m 
tv3m=integ (tv3mdt , tv3mic) 

pv3mdt’ = r/ (v3m*144 . ) * (nv3m*tv3mdt+tv3m*nv3mdt) 
pv3m = integ (pv3mdt , pv3mic) 

ev3mdt=cp3m* ( ( f 3mix+f bypas) *tin- (f 4slup+f 4s2up) *tv3m) 

ii it 

tt dark 4th stage analysis " 

ii it 

f 4slup=rsw(psuc4-p4sl . le . 0 . , 0 . , sqrt (psuc4-p4sl) /k4slup) 

f 4s2up=rsw(psuc4-p4s2 . le . 0 . , 0 . , sqrt (psuc4-p4s2 ) /k4s2up) 

n4sldt = f 4slup-f 4sldn 

n4s2dt = f 4s2up-f4s2dn 

n4sl = p4sl*144 . *v4sl/ (r*t4sl) 

n4s2 = p4s2*144 . *v4s2/ (r*t4s2) 

t4sldt=(e4sldt/cv3m-n4sldt*t4sl) /n4sl 

t4s2dt=(e4s2dt/cv3m-n4s2dt*t4s2 ) /n4s2 

t4sl = integ (t4sldt, t4slic) 

t4s2 = integ (t4s2dt , t4s2ic) 

p4sldt=r/v4sl* (n4sl*t4sldt+t4sl*n4sldt-n4sl*t4sl/v4sl*v4sldt) /144 . 
p4s2dt=r/v4s2 * (n4s2*t4s2dt+t4s2*n4s2dt-n4s2*t4s2 /v4s2 *v4s2dt) /144 . 
p4sl = integ(p4sldt,p4slic) 
p4s2 = integ(p4s2dt,p4s2ic) 

ii ii 

v4sl=(cls4*vef 4-strk/2 . *sin (omega* t) ) *a4sl 
v4s2=(cls4*vef4-strk/2 . *sin (omega *t+pi) ) *a4s2 

ii ii 


v4sldt=-a4sl*strk/2 . *omega*cos (omega*t) 
v4s2dt=-a4s2*strk/2 . *omega*cos (omega*t+pi) 

M It 

e4sldt=cp3m* (f 4slup*tin-f 4sldn*t4sl) -p4sl*v4sldt* 14 4 . / j 
e4s2dt=cp3m* (f 4s2up*tin-f 4s2dn*t4s2 ) -p4s2*v4s2dt*144 . / j 

ti ii 

t4=rsw(f 4sldn.gt . 0. , t4sl-460. , t4) 

H it 

M volume 4 analysis " 

n i» 

f 4sldn=rsw (p4sl-pv4 . le, 0. , 0. , sqrt (p4sl-pv4 ) /k4sldn) 
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f 4s2dn=rsw(p4s2-pv4 • le. 0. ,0. , sqrt (p4s2-pv4 ) /k4s2dn) 

t« " 

nv4dt = (f 4sldn+f 4s2dn) - (f 5slup+f 5s2up) - fcond4 

nv4 = pv4*144 . *v4/ (r*tv4) 

nhg4dt = f4in*x0 - f40ut*x4 

nhg4 = integ(nhg4dt , nhg4ic) 

x4 = nhg4/nv4 

tv4dt = (ev4dt/cv3m-nv4dt*tv4) /nv4 
tv4 = integ(tv4dt , tv4ic) 

pv4dt = r/ (v4*144 • )'* (nv4*tv4dt+tv4*nv4dt) 
pv4 = integ(pv4dt,pv4ic) 

it " 

ev4dt=cp3m* (f 4sldn*t4sl+f 4s2dn*t4s2- (f 5slup+f 5s2up+f cond4 ) *tv4) +qburn4 
tgo4 = rsw(t4.gt.500.0.or.tgo4.eq. 1.0,1. 0,0.0) 
qburn4=hcomb* (f 4sldn+f 4s2dn+. 000116*pv4) *tgo4 
ratio4=(pv4-cp4) / (fp4-cp4) 

f ull4=rsw(pv4-cp4 . le.O. , 0. , fp4*a4/sqrt ( tv4) /cnst) 
f low4=rsw(ratio4 . It, 1 . 0 ,ratio4*full4 , pv4/fp4*full4 ) 


ti »• 

*• 4th stage condensate flow and concentration " 

ti »» 


x4calc=vpl00/pv4 

x4=rsw(pv4 .gt. vplOO. and. x4 calc. It . x3m, x4calc, x3m) 

cp4 = x4*cphg +(l.-x4)*cpa 

cv4 = x4*cvhg +(l.-x4)*cva 

m4 = x4*mhg +(l.-x4)*ma 

f4air = (f 4sldn+f 4s2dn) * ( 1 . -x3m) 

f cond4=f 4air* (x3m* ( 1 . -x4 ) -x4* ( 1 . -x3m) ) / (l.-x3m) / (1. -x4) 

n »' 

" dark 5th stage analysis " 

n »• 

f5slup=rsw(psuc5-p5sl . le. 0 . , 0 . , sqrt (psuc5-p5sl) /k5slup) 

f 5s2up=rsw(psuc5-p5s2 . le. 0. ,0. , sqrt (psuc5-p5s2) /k5s2up) 

nSsldt = f 5slup-f 5sldn 

n5s2dt = f 5s2up-f 5s2dn 

n5sl = p5sl*144.*v5sl/(r*t5sl) 

n5s2 = p5s2*144 . *v5s2/ ( r*t5s2 ) 

t5sldt=(e5sldt/cv4-n5sldt*t5sl) /n5sl 

t5s2dt= (e5s2dt/cv4-n5s2dt*t5s2 ) /n5s2 

t5sl = integ(t5sldt , tSslic) 

t5s2 = integ(t5s2dt , t5s2ic) 

p5sldt=r/v5sl* (n5sl*t5sldt+t5sl*n5sldt-n5sl*t5sl/v5sl*v5sldt) /144 . 
p5s2dt=r/v5s2* (n5s2*t5s2dt+t5s2*n5s2dt-n5s2*t5s2/v5s2*v5s2dt) /144 . 
p5sl = integ(p5sldt f pSslic) 
p5s2 = Integ(p5s2dt,p5s2ic) 

ti » 

v5sl= (cls5*vef 5-strk/2 . *sin (omega *t+pi) ) *a5sl 
v5s2=(cls5*vef 5-strk/2 . *sin (omega*t) ) *a5s2 

it •» 

v5sldt=-a5sl*strk/2 . *omega*cos (omega*t+pi) 
v5s2dt=-a5s2*strk/2 . *omega*cos (omega*t) 

ti •» 

e5sldt=cp4* (f 5slup*tin-f 5sldn*t5sl) -p5sl*v5sldt* 144 . / j 
e5s2dt=cp4* (f 5s2up*tin-f 5s2dn*t5s2) -p5s2*v5s2dt*144 . / j 

it •' 

t5=rsw(f5sldn.gt. 0. , t5sl-460. ,t5) 

II « 

*» volume 5 analysis " 

it " 

f 5sldn=rsw(p5sl-pv5 . le.O. ,0. , sqrt (p5sl-pv5) /k5sldn) 
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f 5s2dn=rsw(p5s2-pv5. le.O. , O, , sqrt (p5s2-pv5) /k5s2dn) 

it 

nv5dt = (f 5sldn+f 5s2dn) - (fvlv+f 5to0+fbypas+fcond5) 

nv5 = pv5*l44.*v5/(r*tv5) 

nhg5dt = f5in*x4 - f5out*x5 

nhg5 = integ(nhg5dt,nhg5ic) 

x5 = nhg5/nv5 

tv5dt = (ev5dt/cv4-nv5dt*tv5) /nv5 
tv5 = integ(tv5dt , tv5ic) 

pv5dt = r/ (v5*144 .) * (nv5*tv5dt+tv5*nv5dt) 
pv5 = integ(pv5dt, pv5ic) 

tt 

ev5dt=cp4* (f5sldn*t5sl+f5s2dn*t5s2- (fvlv+f5to0+fbypas+fcond5) *tv5) . . . 
+qburn5 

tgo5 = rsw(t5 .gt . 500 . 0. or . tgo5 . eq. 1 . 0, 1. 0, 0. 0) 
qburn5=hcomb* (f5sldn+f5s2dn+. 00006*pv5) *tgo5 
ratio5=(pv5-cp5) / (fp5-cp5) 

f ull5=rsw (pv5-cp5 . le . 0 . , 0 . , f p5*a5/sqrt ( tv5) /cnst ) 
f low5=rsw( ratios. It. 1. 0,ratio5*full5,pv5/fp5*full5) 

It 

" volume 5 discharge conditions 

tt 

x5calc = vpl00/pv5 

x5=rsw (pv5 . gt . vplOO . and . x5calc . It . x4 , xScalc , x4 ) 
f5air=(f5sldn+f 5s2dn) * (l.-x4) 

fcond5=f 5air* (x4* (1. -x5) -x5* (1. -x4) ) / (1. -x4) / (1. -x5) 
cp5 = x5*cphg + (l.-x5)*cpa 

" determine suction pressures 

f 12=realpl (1. , ( (f 2slup+f 2s2up) **2 . *losl2* (tin/pvl) *ml) , 0 . ) 

f 23=realpl (1. , ( (f 3slup+f 3s2up) **2 . *los23* (tin/pv2 ) *m2) , 0. ) 

f 34=realpl (1. , ( (f 4slup+f 4s2up) **2 . *los34* (tin/pv3m) *m3m) , 0. ) 

f 45=realpl (1* , ( ( f 5slup+f 5s2up) **2 . *los45* (tin/pv4 ) *m4 ) ,0.) 

tsl=realpl (1. , tl, tslic) 

ts2=realpl (l. , t2,ts2ic) 

ts3=realpl(l. , t3 , ts3ic) 

ts4=realpl (1. ,t4,ts4ic) 

ts5=realpl(l. ,t5,ts5ic) 

liq3=realpl(l. , fcond3, 0. ) 

liq4=realpl(l. , fcond4, 0. ) 

liq5=realpl ( 1 . , f cond5 , 0 . ) 

con3=x3 

con3m=x3m 

con4=x4 

con5=x5 

prl=pfdl/pf si 

pr2=pfd2/pf s2 

pr3=pfd3/pf s3 

pr4=pfd4/pf s4 

pr5=pfd5/pf s5 

psucl = pvO 

psuc2 = pvl-f 12 

psuc3 = pv2-f 23 

psuc4 = pv3m-f34 

psuc5 = pv4-f 45 

pfsl=realpl(l. ,pvO,pfslic) 

pfs2=realpl (1. , psuc2 , pf s2ic) 

pfs3=realpl (1. ,psuc3 ,pfs3ic) 

pfs4=realpl (1. ,psuc4 ,pfs4ic) 

pfs5=realpl(l. f psuc5 , pf s5ic) 

pf dl=realpl ( 1 . ,pvl,pfdlic) 



pfd2=realpl (1 . ,pv2,pfd2ic) 
pfd3=realpl (1. ,pv3 f pfd3ic) 
pf d4=realpl ( 1 . , pv4 , pf d4 ic) 
pfd5=realpl (1. ,pv5,pfd5ic) 
tfdl=realpl(l. , tvl, tfdlic) -460. 
tf d2=realpl ( 1 . ,tv2,tfd2ic) -460. 
tfd3=realpl(l. , tv3 , tfd3ic) -460. 
tfd4=realpl(l. , tv4 , tfd4ic) -460. 
tf d5=realpl ( 1 . , tv5,, tf d5ic) -460 . 


back pressure control valve 

** valve dynamics and controller, assume 1st order 

ydot = (ycmd - y) * tau 

y = limint (ydot, 0.0, 0.0, 1.0) 

open= . 4881*y-l . 2738 *y**2. 4*1. 7857*y**3 . 


dp=pv5-pcbl 

yex=l . 0-0 . 75*dp/pv5 

f vlv=rsw(dp. le . 0 . , 0 . 0 , cg*yex* . 000227 *sqrt (dp*pv5/m0) *open) 

err = pv5-set 

pcmd=kp*err 

icmd=limint (ki*err ,0. 0,0.0, 1.0) 
ycmd=bound (0 . 0 , 1 . 0 , pcmd+icmd) 

it 

n discharge to suction bypass valve 

» valve dynamics and controller, assume 1st order 

ydotl = (ycmdl - yl) * taul 

yl = limint (ydotl , 0 . 0 , 0.0, 1.0) 

openl=. 4881 *y 1-1. 2738 *y 1**2. +1.78 57 *y 1**3 . 

f 5to0 = (pv5) / vcnst 1 / sqrt ( tv5 ) *openl 

errl = rsw(prl-pr 2 . ge . 0 . ,prl-prset , pr2-prset) 

pcmdl=kpl*errl 

icmdl=limint (kil*errl, 0. 0,0.0, 1.0) 
ycmdl=bound (0 .0,1. 0 ,pcmdl+icmdl) 

it 

termt (t.ge. tstp) 
end $"of derivative" 
end $"of dynamic" 
end $"of program" 



lat stag* 



Figure 1 . Model schematic 


Fin, Tin, Cpin 


Fout, Tout 
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STAGE DISCHARGE PRESSURES, PSIA 




































Figure 17. Fifth stage discharge pressure vs. time. 
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Figure 18. Fifth stage relief valve molar flow vs. time. 
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